!
! Simple Fortan90 program that multiplies a matrix and a vector
! using CUBLAS or CUDA
!
! This example is using cudaMallocHost to allocate the matrices in pinned memory.
! Pinned memory enables fast PCI-e transfer.
!
! The code is using the iso_c_binding Fortran 2003 extension.
! It has been tested with the Intel Fortran compiler v10.1 and g95.
!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Interface to cudaMallocHost and cudaFree
module cuda_alloc
 interface
! cudaMallocHost
 integer (C_INT) function cudaMallocHost(buffer, size)  bind(C,name="cudaMallocHost")
  use iso_c_binding
  implicit none
  type (C_PTR)  :: buffer
  integer (C_SIZE_T), value :: size
 end function cudaMallocHost
! cudaFreeHost
 integer (C_INT) function cudaFreeHost(buffer)  bind(C,name="cudaFreeHost")
  use iso_c_binding
  implicit none
  type (C_PTR), value :: buffer
 end function cudaFreeHost
 end interface
end module cuda_alloc
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program matrix_vec_multiply
use iso_c_binding
use cuda_alloc

implicit none

!!!!!!rows and cols must be set to the same value expected in the CUDA/CUBLAS function
integer:: rows = 4096
integer:: cols = 4096

! Define the floating point kind to be  single_precision
integer, parameter :: fp_kind = kind(0.0)

! The allocation is performed by C function call.
! Define the C pointer as type (C_PTR)
 type(C_PTR) :: cptr_A, cptr_B, cptr_C, cptr_CPU

! Define matrices as pointer.
real (fp_kind), dimension(:,:), pointer ::      A
real (fp_kind), dimension(:), pointer ::      B,C,CPU

real ::      time_start,time_end
real ::      time_start1,time_end1

integer:: res

call gpu_init()

! Allocating memory with cudaMallocHost (this way we allocate PINNED/unswappable memory).
! The Fortan arrays, now defined as pointers, are then associated with the C pointers using the
! new interoperability defined in iso_c_binding

 !allocate(A(m1,m1))
 res = cudaMallocHost ( cptr_A, rows*cols*sizeof(fp_kind) )
 call c_f_pointer ( cptr_A, A, (/ rows, cols /) )

 !allocate(B(m1))
 res = cudaMallocHost ( cptr_B, rows*sizeof(fp_kind) )
 call c_f_pointer ( cptr_B, B, (/ rows /) )

 !allocate(C(m1))
 res = cudaMallocHost ( cptr_C, rows*sizeof(fp_kind) )
 call c_f_pointer ( cptr_C, C, (/ rows /) )

  !allocate(CPU(m1))
 res = cudaMallocHost ( cptr_CPU, rows*sizeof(fp_kind) )
 call c_f_pointer ( cptr_CPU, CPU, (/ rows /) )

 call init_host_data (A, B, CPU)

! Compute the matrix product  computation
 !call cpu_time(time_start)
	call cudafunction_cuda (A, B, C)
 !call cpu_time(time_end)

 !call cpu_time(time_start1)
	call cudafunction_CUBLAS (A, B, CPU)
	!we use the vector called CPU instead of creating a new one just for the CUBLAS call
 !call cpu_time(time_end1)

	call vector_eq_test (C, CPU)

! Clean up the C memory allocated. We need to use the cudaFreeHost function call
  res = cudaFreeHost (cptr_A)
  res = cudaFreeHost (cptr_B)
  res = cudaFreeHost (cptr_C)
  res = cudaFreeHost (cptr_CPU)

end program matrix_vec_multiply

